<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<title>Chaotic systems and Lyapunov exponents</title>
<link rel="stylesheet" href="../../../../../../../doc/src/boostbook.css" type="text/css">
<meta name="generator" content="DocBook XSL Stylesheets V1.79.1">
<link rel="home" href="../../index.html" title="Chapter 1. Boost.Numeric.Odeint">
<link rel="up" href="../tutorial.html" title="Tutorial">
<link rel="prev" href="solar_system.html" title="Solar system">
<link rel="next" href="stiff_systems.html" title="Stiff systems">
</head>
<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
<table cellpadding="2" width="100%"><tr>
<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../logo.jpg"></td>
<td align="center"><a href="../../../../../../../index.html">Home</a></td>
<td align="center"><a href="../../../../../../../libs/libraries.htm">Libraries</a></td>
<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
<td align="center"><a href="../../../../../../../more/index.htm">More</a></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="solar_system.html"><img src="../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../tutorial.html"><img src="../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="stiff_systems.html"><img src="../../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
<div class="section">
<div class="titlepage"><div><div><h3 class="title">
<a name="boost_numeric_odeint.tutorial.chaotic_systems_and_lyapunov_exponents"></a><a class="link" href="chaotic_systems_and_lyapunov_exponents.html" title="Chaotic systems and Lyapunov exponents">Chaotic
      systems and Lyapunov exponents</a>
</h3></div></div></div>
<p>
        In this example we present application of odeint to investigation of the
        properties of chaotic deterministic systems. In mathematical terms chaotic
        refers to an exponential growth of perturbations <span class="emphasis"><em>δ x</em></span>.
        In order to observe this exponential growth one usually solves the equations
        for the tangential dynamics which is again an ordinary differential equation.
        These equations are linear but time dependent and can be obtained via
      </p>
<p>
        <span class="emphasis"><em>d δ x / dt = J(x) δ x</em></span>
      </p>
<p>
        where <span class="emphasis"><em>J</em></span> is the Jacobian of the system under consideration.
        <span class="emphasis"><em>δ x</em></span> can also be interpreted as a perturbation of the original
        system. In principle <span class="emphasis"><em>n</em></span> of these perturbations exist,
        they form a hypercube and evolve in the time. The Lyapunov exponents are
        then defined as logarithmic growth rates of the perturbations. If one Lyapunov
        exponent is larger then zero the nearby trajectories diverge exponentially
        hence they are chaotic. If the largest Lyapunov exponent is zero one is usually
        faced with periodic motion. In the case of a largest Lyapunov exponent smaller
        then zero convergence to a fixed point is expected. More information's about
        Lyapunov exponents and nonlinear dynamical systems can be found in many textbooks,
        see for example: E. Ott "Chaos is Dynamical Systems", Cambridge.
      </p>
<p>
        To calculate the Lyapunov exponents numerically one usually solves the equations
        of motion for <span class="emphasis"><em>n</em></span> perturbations and orthonormalizes them
        every <span class="emphasis"><em>k</em></span> steps. The Lyapunov exponent is the average
        of the logarithm of the stretching factor of each perturbation.
      </p>
<p>
        To demonstrate how one can use odeint to determine the Lyapunov exponents
        we choose the Lorenz system. It is one of the most studied dynamical systems
        in the nonlinear dynamics community. For the standard parameters it possesses
        a strange attractor with non-integer dimension. The Lyapunov exponents take
        values of approximately 0.9, 0 and -12.
      </p>
<p>
        The implementation of the Lorenz system is
      </p>
<p>
</p>
<pre class="programlisting"><span class="keyword">const</span> <span class="keyword">double</span> <span class="identifier">sigma</span> <span class="special">=</span> <span class="number">10.0</span><span class="special">;</span>
<span class="keyword">const</span> <span class="keyword">double</span> <span class="identifier">R</span> <span class="special">=</span> <span class="number">28.0</span><span class="special">;</span>
<span class="keyword">const</span> <span class="keyword">double</span> <span class="identifier">b</span> <span class="special">=</span> <span class="number">8.0</span> <span class="special">/</span> <span class="number">3.0</span><span class="special">;</span>

<span class="keyword">typedef</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">array</span><span class="special">&lt;</span> <span class="keyword">double</span> <span class="special">,</span> <span class="number">3</span> <span class="special">&gt;</span> <span class="identifier">lorenz_state_type</span><span class="special">;</span>

<span class="keyword">void</span> <span class="identifier">lorenz</span><span class="special">(</span> <span class="keyword">const</span> <span class="identifier">lorenz_state_type</span> <span class="special">&amp;</span><span class="identifier">x</span> <span class="special">,</span> <span class="identifier">lorenz_state_type</span> <span class="special">&amp;</span><span class="identifier">dxdt</span> <span class="special">,</span> <span class="keyword">double</span> <span class="identifier">t</span> <span class="special">)</span>
<span class="special">{</span>
    <span class="identifier">dxdt</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">=</span> <span class="identifier">sigma</span> <span class="special">*</span> <span class="special">(</span> <span class="identifier">x</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">);</span>
    <span class="identifier">dxdt</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">=</span> <span class="identifier">R</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">x</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span><span class="number">2</span><span class="special">];</span>
    <span class="identifier">dxdt</span><span class="special">[</span><span class="number">2</span><span class="special">]</span> <span class="special">=</span> <span class="special">-</span><span class="identifier">b</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span><span class="number">2</span><span class="special">]</span> <span class="special">+</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span><span class="number">1</span><span class="special">];</span>
<span class="special">}</span>
</pre>
<p>
        We need also to integrate the set of the perturbations. This is done in parallel
        to the original system, hence within one system function. Of course, we want
        to use the above definition of the Lorenz system, hence the definition of
        the system function including the Lorenz system itself and the perturbation
        could look like:
      </p>
<p>
</p>
<pre class="programlisting"><span class="keyword">const</span> <span class="identifier">size_t</span> <span class="identifier">n</span> <span class="special">=</span> <span class="number">3</span><span class="special">;</span>
<span class="keyword">const</span> <span class="identifier">size_t</span> <span class="identifier">num_of_lyap</span> <span class="special">=</span> <span class="number">3</span><span class="special">;</span>
<span class="keyword">const</span> <span class="identifier">size_t</span> <span class="identifier">N</span> <span class="special">=</span> <span class="identifier">n</span> <span class="special">+</span> <span class="identifier">n</span><span class="special">*</span><span class="identifier">num_of_lyap</span><span class="special">;</span>

<span class="keyword">typedef</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">tr1</span><span class="special">::</span><span class="identifier">array</span><span class="special">&lt;</span> <span class="keyword">double</span> <span class="special">,</span> <span class="identifier">N</span> <span class="special">&gt;</span> <span class="identifier">state_type</span><span class="special">;</span>
<span class="keyword">typedef</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">tr1</span><span class="special">::</span><span class="identifier">array</span><span class="special">&lt;</span> <span class="keyword">double</span> <span class="special">,</span> <span class="identifier">num_of_lyap</span> <span class="special">&gt;</span> <span class="identifier">lyap_type</span><span class="special">;</span>

<span class="keyword">void</span> <span class="identifier">lorenz_with_lyap</span><span class="special">(</span> <span class="keyword">const</span> <span class="identifier">state_type</span> <span class="special">&amp;</span><span class="identifier">x</span> <span class="special">,</span> <span class="identifier">state_type</span> <span class="special">&amp;</span><span class="identifier">dxdt</span> <span class="special">,</span> <span class="keyword">double</span> <span class="identifier">t</span> <span class="special">)</span>
<span class="special">{</span>
    <span class="identifier">lorenz</span><span class="special">(</span> <span class="identifier">x</span> <span class="special">,</span> <span class="identifier">dxdt</span> <span class="special">,</span> <span class="identifier">t</span> <span class="special">);</span>

    <span class="keyword">for</span><span class="special">(</span> <span class="identifier">size_t</span> <span class="identifier">l</span><span class="special">=</span><span class="number">0</span> <span class="special">;</span> <span class="identifier">l</span><span class="special">&lt;</span><span class="identifier">num_of_lyap</span> <span class="special">;</span> <span class="special">++</span><span class="identifier">l</span> <span class="special">)</span>
    <span class="special">{</span>
        <span class="keyword">const</span> <span class="keyword">double</span> <span class="special">*</span><span class="identifier">pert</span> <span class="special">=</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">+</span> <span class="number">3</span> <span class="special">+</span> <span class="identifier">l</span> <span class="special">*</span> <span class="number">3</span><span class="special">;</span>
        <span class="keyword">double</span> <span class="special">*</span><span class="identifier">dpert</span> <span class="special">=</span> <span class="identifier">dxdt</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">+</span> <span class="number">3</span> <span class="special">+</span> <span class="identifier">l</span> <span class="special">*</span> <span class="number">3</span><span class="special">;</span>
        <span class="identifier">dpert</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">=</span> <span class="special">-</span> <span class="identifier">sigma</span> <span class="special">*</span> <span class="identifier">pert</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">+</span> <span class="number">10.0</span> <span class="special">*</span> <span class="identifier">pert</span><span class="special">[</span><span class="number">1</span><span class="special">];</span>
        <span class="identifier">dpert</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">=</span> <span class="special">(</span> <span class="identifier">R</span> <span class="special">-</span> <span class="identifier">x</span><span class="special">[</span><span class="number">2</span><span class="special">]</span> <span class="special">)</span> <span class="special">*</span> <span class="identifier">pert</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">pert</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">*</span> <span class="identifier">pert</span><span class="special">[</span><span class="number">2</span><span class="special">];</span>
        <span class="identifier">dpert</span><span class="special">[</span><span class="number">2</span><span class="special">]</span> <span class="special">=</span> <span class="identifier">x</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">*</span> <span class="identifier">pert</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">+</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">*</span> <span class="identifier">pert</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">b</span> <span class="special">*</span> <span class="identifier">pert</span><span class="special">[</span><span class="number">2</span><span class="special">];</span>
    <span class="special">}</span>
<span class="special">}</span>
</pre>
<p>
      </p>
<p>
        The perturbations are stored linearly in the <code class="computeroutput"><span class="identifier">state_type</span></code>
        behind the state of the Lorenz system. The problem of lorenz() and lorenz_with_lyap() having different
        state types may be solved putting the Lorenz system inside a functor with
        templatized arguments:
      </p>
<p>
</p>
<pre class="programlisting"><span class="keyword">struct</span> <span class="identifier">lorenz</span>
<span class="special">{</span>
    <span class="keyword">template</span><span class="special">&lt;</span> <span class="keyword">class</span> <span class="identifier">StateIn</span> <span class="special">,</span> <span class="keyword">class</span> <span class="identifier">StateOut</span> <span class="special">,</span> <span class="keyword">class</span> <span class="identifier">Value</span> <span class="special">&gt;</span>
    <span class="keyword">void</span> <span class="keyword">operator</span><span class="special">()(</span> <span class="keyword">const</span> <span class="identifier">StateIn</span> <span class="special">&amp;</span><span class="identifier">x</span> <span class="special">,</span> <span class="identifier">StateOut</span> <span class="special">&amp;</span><span class="identifier">dxdt</span> <span class="special">,</span> <span class="identifier">Value</span> <span class="identifier">t</span> <span class="special">)</span>
    <span class="special">{</span>
        <span class="identifier">dxdt</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">=</span> <span class="identifier">sigma</span> <span class="special">*</span> <span class="special">(</span> <span class="identifier">x</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">);</span>
        <span class="identifier">dxdt</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">=</span> <span class="identifier">R</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">x</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span><span class="number">2</span><span class="special">];</span>
        <span class="identifier">dxdt</span><span class="special">[</span><span class="number">2</span><span class="special">]</span> <span class="special">=</span> <span class="special">-</span><span class="identifier">b</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span><span class="number">2</span><span class="special">]</span> <span class="special">+</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span><span class="number">1</span><span class="special">];</span>
    <span class="special">}</span>
<span class="special">};</span>

<span class="keyword">void</span> <span class="identifier">lorenz_with_lyap</span><span class="special">(</span> <span class="keyword">const</span> <span class="identifier">state_type</span> <span class="special">&amp;</span><span class="identifier">x</span> <span class="special">,</span> <span class="identifier">state_type</span> <span class="special">&amp;</span><span class="identifier">dxdt</span> <span class="special">,</span> <span class="keyword">double</span> <span class="identifier">t</span> <span class="special">)</span>
<span class="special">{</span>
    <span class="identifier">lorenz</span><span class="special">()(</span> <span class="identifier">x</span> <span class="special">,</span> <span class="identifier">dxdt</span> <span class="special">,</span> <span class="identifier">t</span> <span class="special">);</span>
    <span class="special">...</span>
<span class="special">}</span>
</pre>
<p>
        This works fine and <code class="computeroutput"><span class="identifier">lorenz_with_lyap</span></code>
        can be used for example via
</p>
<pre class="programlisting"><span class="identifier">state_type</span> <span class="identifier">x</span><span class="special">;</span>
<span class="comment">// initialize x..</span>

<span class="identifier">explicit_rk4</span><span class="special">&lt;</span> <span class="identifier">state_type</span> <span class="special">&gt;</span> <span class="identifier">rk4</span><span class="special">;</span>
<span class="identifier">integrate_n_steps</span><span class="special">(</span> <span class="identifier">rk4</span> <span class="special">,</span> <span class="identifier">lorenz_with_lyap</span> <span class="special">,</span> <span class="identifier">x</span> <span class="special">,</span> <span class="number">0.0</span> <span class="special">,</span> <span class="number">0.01</span> <span class="special">,</span> <span class="number">1000</span> <span class="special">);</span>
</pre>
<p>
        This code snippet performs 1000 steps with constant step size 0.01.
      </p>
<p>
        A real world use case for the calculation of the Lyapunov exponents of Lorenz
        system would always include some transient steps, just to ensure that the
        current state lies on the attractor, hence it would look like
      </p>
<p>
</p>
<pre class="programlisting"><span class="identifier">state_type</span> <span class="identifier">x</span><span class="special">;</span>
<span class="comment">// initialize x</span>
<span class="identifier">explicit_rk4</span><span class="special">&lt;</span> <span class="identifier">state_type</span> <span class="special">&gt;</span> <span class="identifier">rk4</span><span class="special">;</span>
<span class="identifier">integrate_n_steps</span><span class="special">(</span> <span class="identifier">rk4</span> <span class="special">,</span> <span class="identifier">lorenz</span> <span class="special">,</span> <span class="identifier">x</span> <span class="special">,</span> <span class="number">0.0</span> <span class="special">,</span> <span class="number">0.01</span> <span class="special">,</span> <span class="number">1000</span> <span class="special">);</span>
</pre>
<p>
        The problem is now, that <code class="computeroutput"><span class="identifier">x</span></code>
        is the full state containing also the perturbations and <code class="computeroutput"><span class="identifier">integrate_n_steps</span></code>
        does not know that it should only use 3 elements. In detail, odeint and its
        steppers determine the length of the system under consideration by determining
        the length of the state. In the classical solvers, e.g. from Numerical Recipes,
        the problem was solved by pointer to the state and an appropriate length,
        something similar to
      </p>
<p>
</p>
<pre class="programlisting"><span class="keyword">void</span> <span class="identifier">lorenz</span><span class="special">(</span> <span class="keyword">double</span><span class="special">*</span> <span class="identifier">x</span> <span class="special">,</span> <span class="keyword">double</span> <span class="special">*</span><span class="identifier">dxdt</span> <span class="special">,</span> <span class="keyword">double</span> <span class="identifier">t</span><span class="special">,</span> <span class="keyword">void</span><span class="special">*</span> <span class="identifier">params</span> <span class="special">)</span>
<span class="special">{</span>
    <span class="special">...</span>
<span class="special">}</span>

<span class="keyword">int</span> <span class="identifier">system_length</span> <span class="special">=</span> <span class="number">3</span><span class="special">;</span>
<span class="identifier">rk4</span><span class="special">(</span> <span class="identifier">x</span> <span class="special">,</span> <span class="identifier">system_length</span> <span class="special">,</span> <span class="identifier">t</span> <span class="special">,</span> <span class="identifier">dt</span> <span class="special">,</span> <span class="identifier">lorenz</span> <span class="special">);</span>
</pre>
<p>
      </p>
<p>
        But odeint supports a similar and much more sophisticated concept: <a href="http://www.boost.org/doc/libs/release/libs/range/" target="_top">Boost.Range</a>.
        To make the steppers and the system ready to work with <a href="http://www.boost.org/doc/libs/release/libs/range/" target="_top">Boost.Range</a>
        the system has to be changed:
      </p>
<p>
</p>
<pre class="programlisting"><span class="keyword">struct</span> <span class="identifier">lorenz</span>
<span class="special">{</span>
    <span class="keyword">template</span><span class="special">&lt;</span> <span class="keyword">class</span> <span class="identifier">State</span> <span class="special">,</span> <span class="keyword">class</span> <span class="identifier">Deriv</span> <span class="special">&gt;</span>
    <span class="keyword">void</span> <span class="keyword">operator</span><span class="special">()(</span> <span class="keyword">const</span> <span class="identifier">State</span> <span class="special">&amp;</span><span class="identifier">x_</span> <span class="special">,</span> <span class="identifier">Deriv</span> <span class="special">&amp;</span><span class="identifier">dxdt_</span> <span class="special">,</span> <span class="keyword">double</span> <span class="identifier">t</span> <span class="special">)</span> <span class="keyword">const</span>
    <span class="special">{</span>
        <span class="keyword">typename</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">range_iterator</span><span class="special">&lt;</span> <span class="keyword">const</span> <span class="identifier">State</span> <span class="special">&gt;::</span><span class="identifier">type</span> <span class="identifier">x</span> <span class="special">=</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">begin</span><span class="special">(</span> <span class="identifier">x_</span> <span class="special">);</span>
        <span class="keyword">typename</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">range_iterator</span><span class="special">&lt;</span> <span class="identifier">Deriv</span> <span class="special">&gt;::</span><span class="identifier">type</span> <span class="identifier">dxdt</span> <span class="special">=</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">begin</span><span class="special">(</span> <span class="identifier">dxdt_</span> <span class="special">);</span>

        <span class="identifier">dxdt</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">=</span> <span class="identifier">sigma</span> <span class="special">*</span> <span class="special">(</span> <span class="identifier">x</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">);</span>
        <span class="identifier">dxdt</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">=</span> <span class="identifier">R</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">x</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span><span class="number">2</span><span class="special">];</span>
        <span class="identifier">dxdt</span><span class="special">[</span><span class="number">2</span><span class="special">]</span> <span class="special">=</span> <span class="special">-</span><span class="identifier">b</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span><span class="number">2</span><span class="special">]</span> <span class="special">+</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span><span class="number">1</span><span class="special">];</span>
    <span class="special">}</span>
<span class="special">};</span>
</pre>
<p>
      </p>
<p>
        This is in principle all. Now, we only have to call <code class="computeroutput"><span class="identifier">integrate_n_steps</span></code>
        with a range including only the first 3 components of <span class="emphasis"><em>x</em></span>:
      </p>
<p>
</p>
<pre class="programlisting"><span class="comment">// explicitly choose range_algebra to override default choice of array_algebra</span>
<span class="identifier">runge_kutta4</span><span class="special">&lt;</span> <span class="identifier">state_type</span> <span class="special">,</span> <span class="keyword">double</span> <span class="special">,</span> <span class="identifier">state_type</span> <span class="special">,</span> <span class="keyword">double</span> <span class="special">,</span> <span class="identifier">range_algebra</span> <span class="special">&gt;</span> <span class="identifier">rk4</span><span class="special">;</span>

<span class="comment">// perform 10000 transient steps</span>
<span class="identifier">integrate_n_steps</span><span class="special">(</span> <span class="identifier">rk4</span> <span class="special">,</span> <span class="identifier">lorenz</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">make_pair</span><span class="special">(</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">+</span> <span class="identifier">n</span> <span class="special">)</span> <span class="special">,</span> <span class="number">0.0</span> <span class="special">,</span> <span class="identifier">dt</span> <span class="special">,</span> <span class="number">10000</span> <span class="special">);</span>
</pre>
<p>
      </p>
<div class="note"><table border="0" summary="Note">
<tr>
<td rowspan="2" align="center" valign="top" width="25"><img alt="[Note]" src="../../../../../../../doc/src/images/note.png"></td>
<th align="left">Note</th>
</tr>
<tr><td align="left" valign="top"><p>
          Note that when using <a href="http://www.boost.org/doc/libs/release/libs/range/" target="_top">Boost.Range</a>,
          we have to explicitly configure the stepper to use the <code class="computeroutput"><span class="identifier">range_algebra</span></code>
          as otherwise odeint would automatically chose the <code class="computeroutput"><span class="identifier">array_algebra</span></code>,
          which is incompatible with the usage of <a href="http://www.boost.org/doc/libs/release/libs/range/" target="_top">Boost.Range</a>,
          because the original state_type is an <code class="computeroutput"><span class="identifier">array</span></code>.
        </p></td></tr>
</table></div>
<p>
        Having integrated a sufficient number of transients steps we are now able
        to calculate the Lyapunov exponents:
      </p>
<div class="orderedlist"><ol class="orderedlist" type="1">
<li class="listitem">
            Initialize the perturbations. They are stored linearly behind the state
            of the Lorenz system. The perturbations are initialized such that <span class="emphasis"><em>p
            <sub>​ij</sub> = δ <sub>​ij</sub></em></span>, where <span class="emphasis"><em>p <sub>​ij</sub></em></span> is the <span class="emphasis"><em>j</em></span>-component
            of the <span class="emphasis"><em>i</em></span>.-th perturbation and <span class="emphasis"><em>δ <sub>​ij</sub></em></span>
            is the Kronecker symbol.
          </li>
<li class="listitem">
            Integrate 100 steps of the full system with perturbations
          </li>
<li class="listitem">
            Orthonormalize the perturbation using Gram-Schmidt orthonormalization
            algorithm.
          </li>
<li class="listitem">
            Repeat step 2 and 3. Every 10000 steps write the current Lyapunov exponent.
          </li>
</ol></div>
<p>
</p>
<pre class="programlisting"><span class="identifier">fill</span><span class="special">(</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()+</span><span class="identifier">n</span> <span class="special">,</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">end</span><span class="special">()</span> <span class="special">,</span> <span class="number">0.0</span> <span class="special">);</span>
<span class="keyword">for</span><span class="special">(</span> <span class="identifier">size_t</span> <span class="identifier">i</span><span class="special">=</span><span class="number">0</span> <span class="special">;</span> <span class="identifier">i</span><span class="special">&lt;</span><span class="identifier">num_of_lyap</span> <span class="special">;</span> <span class="special">++</span><span class="identifier">i</span> <span class="special">)</span> <span class="identifier">x</span><span class="special">[</span><span class="identifier">n</span><span class="special">+</span><span class="identifier">n</span><span class="special">*</span><span class="identifier">i</span><span class="special">+</span><span class="identifier">i</span><span class="special">]</span> <span class="special">=</span> <span class="number">1.0</span><span class="special">;</span>
<span class="identifier">fill</span><span class="special">(</span> <span class="identifier">lyap</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">lyap</span><span class="special">.</span><span class="identifier">end</span><span class="special">()</span> <span class="special">,</span> <span class="number">0.0</span> <span class="special">);</span>

<span class="keyword">double</span> <span class="identifier">t</span> <span class="special">=</span> <span class="number">0.0</span><span class="special">;</span>
<span class="identifier">size_t</span> <span class="identifier">count</span> <span class="special">=</span> <span class="number">0</span><span class="special">;</span>
<span class="keyword">while</span><span class="special">(</span> <span class="keyword">true</span> <span class="special">)</span>
<span class="special">{</span>

    <span class="identifier">t</span> <span class="special">=</span> <span class="identifier">integrate_n_steps</span><span class="special">(</span> <span class="identifier">rk4</span> <span class="special">,</span> <span class="identifier">lorenz_with_lyap</span> <span class="special">,</span> <span class="identifier">x</span> <span class="special">,</span> <span class="identifier">t</span> <span class="special">,</span> <span class="identifier">dt</span> <span class="special">,</span> <span class="number">100</span> <span class="special">);</span>
    <span class="identifier">gram_schmidt</span><span class="special">&lt;</span> <span class="identifier">num_of_lyap</span> <span class="special">&gt;(</span> <span class="identifier">x</span> <span class="special">,</span> <span class="identifier">lyap</span> <span class="special">,</span> <span class="identifier">n</span> <span class="special">);</span>
    <span class="special">++</span><span class="identifier">count</span><span class="special">;</span>

    <span class="keyword">if</span><span class="special">(</span> <span class="special">!(</span><span class="identifier">count</span> <span class="special">%</span> <span class="number">100000</span><span class="special">)</span> <span class="special">)</span>
    <span class="special">{</span>
        <span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="identifier">t</span><span class="special">;</span>
        <span class="keyword">for</span><span class="special">(</span> <span class="identifier">size_t</span> <span class="identifier">i</span><span class="special">=</span><span class="number">0</span> <span class="special">;</span> <span class="identifier">i</span><span class="special">&lt;</span><span class="identifier">num_of_lyap</span> <span class="special">;</span> <span class="special">++</span><span class="identifier">i</span> <span class="special">)</span> <span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"\t"</span> <span class="special">&lt;&lt;</span> <span class="identifier">lyap</span><span class="special">[</span><span class="identifier">i</span><span class="special">]</span> <span class="special">/</span> <span class="identifier">t</span> <span class="special">;</span>
        <span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
    <span class="special">}</span>
<span class="special">}</span>
</pre>
<p>
      </p>
<p>
        The full code can be found here: <a href="https://github.com/headmyshoulder/odeint-v2/blob/master/examples/chaotic_system.cpp" target="_top">chaotic_system.cpp</a>
      </p>
</div>
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
<td align="left"></td>
<td align="right"><div class="copyright-footer">Copyright © 2009-2015 Karsten Ahnert and Mario Mulansky<p>
        Distributed under the Boost Software License, Version 1.0. (See accompanying
        file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>)
      </p>
</div></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="solar_system.html"><img src="../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../tutorial.html"><img src="../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="stiff_systems.html"><img src="../../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
</body>
</html>
